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Abstract 

We consider the problem of optimal reactive power compensation for the minimization of power distribution 
losses in a smart microgrid. We first propose an approximate model for the power distribution network, which 
allows us to cast the problem into the class of convex quadratic, linearly constrained, optimization problems. Then, 
we design a randomized, gossip-like optimization algorithm based on that model. We show how a distributed 
approach is possible, where agents have a partial knowledge of the problem parameters and state, and can only 
perform local measurements. For the proposed algorithm, we provide conditions for convergence together with an 
analytic characterization of the convergence speed. The analysis shows that the best performance can be achieved 
when we command cooperation among agents that are neighbors in the electric topology. Numerical simulations 
are included to validate the proposed model and to confirm the analytic results about the performance of the 
proposed algorithm. 
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I. Introduction 

Most of the distributed optimization methods have been derived for the problem of dispatching part of a large 
scale optimization algorithm to different processing units [1]. When the same methods are applied to networked 
control systems (NCS) [2], however, different issues arise: the way in which decision variables are assigned 
to different agents is not part of the designer degrees of freedom; agents have a local and limited knowledge 
of the problem parameters and of the system state; moreover, the information exchange between agents can 
occur not only via a given communication channel, but also via local actuation and subsequent measurement 
performed on an underlying physical system. 

The extent of these issues depends on the particular application. In this work we present a specific scenario, 
belonging to the motivating framework of smart electrical power distribution networks [3], [4], in which these 
features play a central role. 

In the last decade, the introduction of distributed microgeneration of electric energy (driven by technological 
advances, economical reasons, and environmental issues), together with an increased demand and the need for 
higher quality of service, has been driving the integration of a large amount of information and communication 
technologies (ICT) into the power distribution network. 

Among the many different aspects of this transition, we focus on the control of the microgenerators inside a 
smart microgrid [5], [6]. A microgrid is a portion of the low-voltage power distribution network that is managed 
autonomously from the rest of the network, in order to achieve better quality of the service, improve efficiency, 
and pursue specific economic interests. Together with the loads connected to the microgrid (both residential 
and industrial customers), we also have microgeneration devices (solar panels, combined heat-and-power plants, 
micro wind turbines, etc.). These devices are connected to the microgrid via electronic interfaces (inverters), 
whose main task is to enable the injection of the produced power into the microgrid. 

However, these devices can also perform different other tasks, denoted as ancillary services [7], [8], [9]: 
reactive power compensation, harmonic compensation, voltage support. 

In this work we consider the problem of optimal reactive power compensation. Loads belonging to the 
microgrid may require a sinusoidal current which is not in phase with voltage. A convenient description for 
this, consists in saying that they demand reactive power together with active power, associated with out-of-phase 
and in-phase components of the current, respectively. Reactive power is not a "real" physical power, meaning 
that there is no energy conversion involved nor fuel costs to produce it. Like active power flows, reactive power 
flows contribute to power losses on the lines, cause voltage drop, and may lead to grid instability. It is therefore 
preferable to minimize reactive power flows by producing it as close as possible to the users that need it. 

We explore the possibility of using the electronic interface of the microgeneration units to optimize the flows 
of reactive power in the microgrid. Indeed, the inverters of these units are generally oversized, because most 
of the distributed energy sources are intermittent in time, and the electronic interface is designed according to 
the peak power production. When they are not working at the rated power, these inverters can be commanded 
to inject a desired amount of reactive power at no cost [10]. 

This idea has been recently investigated in the literature on power systems [11], [12], [13]. However, these 
works consider a centralized scenario in which the parameters of the entire power grid are known, the controller 
can access the entire system state, and the few agents deployed in the distribution network are dispatched by a 
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central processing unit. The contribution of this work consists in casting the same problem in the framework of 
networked control systems, distributed optimization, and large scale complex systems. This approach allows the 
design of algorithms and solutions which can guarantee scalability, robustness to node insertion and removal, and 
compliance with the actual communication capabilities of the agents. Up to now, the few attempts of applying 
these tools to the power distribution networks have focused on grids comprising a large number of mechanical 
synchronous generators, instead of power inverters. See for example the stability analysis for these systems 
in [14] and the decentralized control synthesis in [15]. Instead, preliminary attempts of designing distributed 
reactive power compensation strategies for power inverters have been performed only very recently. However, 
the available results in this sense are mainly supported by simulations and heuristic approaches [16], and in 
some cases do not implement any communication or synergistic behavior of the agents [17]. 

The contribution of this paper is twofold. On one side, it consists in proposing a rigorous analytic derivation 
of an approximate model of the power flows which is commonly adopted in the power system literature and 
which allows to cast the optimal reactive power flow (ORPF) problem into a quadratic optimization (Section 
III). The second contribution is to propose and to analyze a distributed optimization algorithm whose operation 
is strongly based on the developed model (Section IV). In Section V we show under which conditions the 
previous algorithm converges to the global optimal solution. In Section VI we show how this algorithm can be 
optimized by a proper choice of the communication strategy. In Section VII, we finally validate the proposed 
model and the proposed algorithm via simulations. 

II. Mathematical preliminaries and notation 

Let Q — (V, £, cr, t) be a directed graph, where V is the set of nodes, £ is the set of edges, and cr, r : £ — > V 
are two functions such that edge e e £ goes from the source node er(e) to the terminal node r(e). Two edges 
e and e' are consecutive if the intersection {a(e),r(e)} n {cr(e'), -r(e')} is not empty. A path is a sequence of 
consecutive edges (regardless of their direction). 

In the rest of the paper we will often introduce complex-valued functions defined on the nodes and on the 
edges. These functions will also be intended as vectors in C" and C r (where n — \V\,r — \£\). We denote by 
" the operation of element-wise complex conjugation, and by T the operation of matrix transposition. 

Let moreover A <E {0, ±l} rx ™ be the incidence matrix of the graph Q, defined via its elements 

— 1 if v = cr(e) 
[A] ev = < 1 if v = r(e) 
otherwise. 

V 

If W is a subset of indices, we define by lyy the column vector whose elements are 



[lw]» < 



1 if v e W 
otherwise. 

Similarly, if w is an index, we denote by l w the column vector whose value is 1 in position w, and elsewhere, 
and we denote by 1 the column vector of all ones. Notice that, if the graph Q is connected (i.e. for every pair 
of nodes v, w e V there is a path connecting them), then 1 is the only vector in ker A. 
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Figure 1. Schematic representation of the adopted microgrid model. In the lower panel a circuit representation is given, where black 
diamonds are microgenerators, white diamonds are loads, and the left-most element of the circuit represents the PCC. The middle panel 
illustrates the adopted graph representation for the same microgrid. Circled nodes represent microgenerators and the edge orientation is 
consistent with the sign of the power line currents. The upper panel shows how the compensators can be divided into overlapping clusters 
in order to implement the control algorithm proposed in Section IV. Each microgenerator is capable of measuring and actuating the system 
locally, while each cluster is also provided with some computational capability. 

An undirected graph Q is a graph in which for every edge e e £, there exists an edge e' € £ such that 
cr(e') = r(e) and r(e') = <r(e). In case the graph has no multiple edges and no self loops, we can also 
describe the edges of an undirected graph as subsets e C V with cardinality |e| = 2. By extension, we define 
a hypergraph % as a pair (V,£) in which edges are subsets of V of arbitrary cardinality. 

III. Problem formulation 

A. A model of a microgrid 

We define a smart microgrid as a portion of the power distribution network that is connected to the power 
transmission network in one point, and hosts a number of loads and micro power generators, as described for 
example in [5], [6]. 
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For the purpose of this paper, we model a microgrid as a directed graph Q, in which edges represent the 
power lines, and nodes represent both loads and generators that are connected to the microgrid. These include 
the residential and industrial consumers, microgenerators, and also the point of connection of the microgrid to 
the transmission grid (called point of common coupling, or PCC). 

We limit our study to the steady state behavior of the system, when all voltages and currents are sinusoidal 
signals at the same frequency cj n . Each signal can therefore be represented via a complex number whose absolute 
value corresponds to the signal root-mean-square value, and whose phase corresponds to the phase of the signal 
at t = 0. Therefore the complex number y — \y\e^ y represents the signal y(t) = \y\^/2 sin(uj t + Zy). 

In this notation, the steady state of a microgrid is described by the following system variables: 

• ii£ C", where u v is the grid voltage at node v; 

• ieC™, where i v is the current injected by node v into the grid; 

• £ G C, where £ e is the current flowing on the edge e. 
The following constraints are satisfied by u, i and £ 

A T £ + i = 0, (1) 
Au + Zi = 0, (2) 

where A is the incidence matrix of Q, and Z = diag(z e ,e e £) is the diagonal matrix of line impedances, 
z e being the impedance of the microgrid power line corresponding to the edge e. Equation (1) corresponds to 
Kirchhoff's current law (KCL) at the nodes, while (2) describes the voltage drop on the edges of the graph. 

Each node v of the microgrid is then characterized by a law relating its injected current i v with its voltage u v . 
We model the node corresponding to the PCC (which we assume to be node 0) as a constant voltage generator, 
i.e. 

u Q = U . (3) 

We assume instead that the voltage u v and the current i v of every node v, but the PCC, satisfy the following 
law 

, V«eV\{0}, (4) 

where s v is the nominal complex power and r] v is a characteristic parameter of the node v. The model (4) is 
called exponential model [18] and is widely adopted in the literature on power flow analysis [19]. Notice that 
s v is the complex power that the node would inject into the grid if the voltage at its point of connection was 
Uq. Notice also that, by properly choosing the parameters r/ v , equation (4) describes the behavior of constant 
power, constant current, and constant impedance devices (with i] v = 0,1,2 respectively). In this sense, this 
model is also a generalization of ZIP models [18], which also are very common in the power system literature. 

It is correct to say that equation (4) models the steady state of the vast majority of residential and industrial 
loads. Microgenerators also fit in this model, with r\ v = 0, as they generally are commanded via a complex 
power reference and they are capable of injecting it independently from the voltage at their point of connection 
[5], [6]. 

The task of solving the system of nonlinear equations given by (1), (2), (3), and (4) to obtain the grid voltages 
and currents, given the network parameters, the injected nominal powers {s v ,v € V\{0}} at every node, and 



u v i v 



U 
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the nominal voltage U at the PCC, has been extensively covered in the literature under the denomination of 
power flow analysis (see for example [20, Chapter 3]). 

In the following, we derive an approximate model for the microgrid state, which will be used later for the 
setup of the optimization problem and for the derivation of the proposed distributed algorithm. To do so, a 
couple of technical lemmas are needed. 

Lemma 1. Let L be the complex valued Laplacian L := A T Z~ l A. There exists a unique symmetric matrix 
X eC nxn such that 

XL = I-\\l 
X1 = 0, 

where, we recall, [1q] v — 1 for v = 0, and otherwise and I is the identity matrix. 



(5) 



Proof: Let us first prove the existence of X. Note that kerX = spanl = ker(7 — 11q), and therefore 
there exists X' E C nxn such that X L = (I - lljf). Let X = X'(I - 1 1 T ). Then 

XL = X'(I - 1 1 T )L = X L = I - 11 T , 

X1 Q = X'(/-l l T )lo = 0. 

Existence is therefore guaranteed. To prove uniqueness, notice that 



Therefore 



and uniqueness of X follows from the uniqueness of the inverse. Moreover, as L = L , we have 



X 


1 




L 


lo 




xl + lis 


Xl 


1 T 







1 T 
1 o 









1 T L 


l T lo 








X 


1 




L 1 


-1 








1 T 







lo 







I 










1 



~X T 


1 




L 


lo 


-T 


L T 1 


-l 


X 


1 


1 T 







■■-o 







IT 0_ 




1 T 






and therefore X = X . ■ 
The matrix X depends only on the topology of the microgrid power lines and on their impedance (compare 
it with the definition of Green matrix in [21]). Indeed, it can be shown that, for every pair of nodes (u,v), 

(1„ - 1 V ) T X(1 U - l v ) = Zf v , (6) 

where Z eS is a complex-valued matrix whose elements Z^l represent the effective impedance of the power 
lines between node u and v. 

If we fix the nominal voltage Uq at the PCC, all the currents i and the voltages u are therefore determined 
by the equations 

u = Xi + U 1 

l T i = (7) 



u v 
Uo 



, Vv e V\{0}, 
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where the first equation results from (1), (2), and (3) together with Lemma 1. 

We can see the currents i and the voltages u as functions i(Uo),u(Uo) of Uq. The following proposition 
provides the Taylor approximation of i(Uo) for large Uo- 

Proposition 2. Let s := — 2~2 v ev\{o} s v Then for all v e V we have that 

iv(U ) = (s v + 5 v (U ))±- (8) 
where S v (Uo) is infinitesimal when Uo tends to infinity. 

The proof of this proposition is based on elementary multivariable analysis but requires a quite involved 
notation. For this reason it is given in the appendix. The quality of the approximation proposed in the previous 
proposition relies on having large nominal voltage Uq and relatively small currents injected by the inverters 
(or supplied to the loads). This assumption corresponds to correct design and operation of power distribution 
networks, where indeed the nominal voltage is chosen sufficiently large (subject to other functional constraints) 
in order to deliver electric power to the loads with relatively small power losses on the power lines. 

Numerical simulation in Section VII will indeed show that the approximation is extremely good when power 
distribution networks are operated in their regular regime. 

Remark. Notice that this approximation, in similar fashions, has been used in the literature before for the 
problem of estimating power flows on the power lines (see among the others [22], [17] and references therein). 
It is also possible to show the analogy between this approximation and the quite common assumption of power 
flow conservation across the power network. The analytical rigorous justification proposed here, however, 
allows us to estimate the approximation error, and, more importantly, will provide the tools to understand what 
information on the system state can be gathered by properly sensing the microgrid. 

B. Power losses minimization problem 

Similarly to what has been done for example in [23], we choose active power losses on the power lines as 
a metric for optimality of reactive power flows. The total active power losses on the edges are given by 

J ' = l^ e )! 2 Rez ( e ) = iT Re(i)u, 

eeS 

where we used (2) to express £ as a function of u. Via (7) and by exploiting the properties (5) of X, we have 

J' = i T X Re(L)Xi = i T Re(X)i. 

Let s be the vector of all nominal complex powers s v , including s := — X^ev\{o} Sv - Then, by plugging in 
the approximate system state (8), we have 

J' = ^ T MX)s + ^J(U ,s) 



1 jP T Re{X)p+ — l -^q T Rc(X)q+ — l —J(U ,s) 



where s is decomposed into the injected active power p = Re(s) and the injected reactive power q = Im(s), 
and where J(Uo, s) := S T Rc(X)s + s T Re(X)S + S T Re(X)S is infinitesimal as Uo tends to infinity and so 
it can be neglected. 
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Therefore, via the approximate model (8), we have been able to 

• approximate power losses as a quadratic function of the injected power; 

• decouple the problem of optimal power flows into the problem of optimal active and reactive power 
injection. 

In the scenario that we are considering, we are allowed to command only a subset C C V of the nodes of the 
microgrid (namely the microgenerators, also named compensators in this framework). We denote by m := \C\ its 
cardinality. Moreover, we assume that for these agents we are only allowed to command the amount of reactive 
power injected into the grid, as the decision on the amount of active power follows imperative economic criteria 
(for example, in the case of renewable energy sources, any available active power is typically injected into the 
grid to replace generation from traditional plants, which are more expensive and exhibit a worse environmental 
impact - see [10]). 

The problem of optimal reactive power injection at the compensators can therefore be expressed as a quadratic, 
linearly constrained problem, in the form 

min J(q), where J(q) = -q T Re(X)q, 
q 2 

subject to l T q = ( 9 ) 

q v = Im(s v ) Vv e V\C, 

{lm(s v ),v e V\C} being the nominal amount of reactive power injected by the nodes that cannot be com- 
manded. 

The solution of the optimization problem (9) would not pose any challenge if the nodes knew the problem 
parameters X and {lm(s v ),v e V\C}. These quantities depend both on the grid parameters and the power 
demand of all the microgrid devices, therefore it is unpractical that every agent is capable of retrieving all this 
information, and, in the case of a large scale system, it would be also preferable to avoid the presence of a 
centralized agent collecting all the necessary data. 

IV. A RANDOMIZED DISTRIBUTED ALGORITHM 

In this section, we present an algorithm that can be distributed across the agents of the microgrid, meaning that 
it can be implemented by the compensators without the supervision of any central controller, by performing local 
measurements, data processing, and actuation. The proposed approach consists in decomposing the optimization 
problem into smaller, tractable subproblems that are assigned to small groups of agents. We will show that the 
agents can solve their corresponding subproblems via local measurement, local knowledge of the grid topology, 
and minimal communication. We will finally show that the iterative, randomized, solution of these subproblems 
yields the solution of the original optimization problem. 

It is worth remarking that the decomposition methods proposed in most of the literature on distributed 
optimization, for example in [24], cannot be applied to this problem because the cost function in (9) is not 
separable into a sum of individual terms for each agent. Approaching the decomposition of this optimization 
problem via its dual formulation (as proposed in many works including the recent [25]) is also unlikely to 
succeed, as the feasibility of the system state must be guaranteed at any time during the optimization process. 
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A. Optimization problem decomposition 

Let the compensators be divided into I possibly overlapping sets C\, . . . , Ce, with Ui=i = TM S family 
of subsets can be interpreted as the edges of a hypergraph W defined over the set of compensator nodes 
C. Nodes belonging to the same set (or cluster) are able to communicate each other, and they are therefore 
capable of coordinating their actions and sharing their measurements. Each cluster is also provided with some 
computational capabilities, in the form of a cluster supervisor, possibly one of the compensators. 

The proposed optimization algorithm consists of the following, repeated steps: 

1) a set Ci{t) is chosen at a certain discrete time t = 0, 1, 2, . . . where i(t) e {1, ...,£}; 

2) the agents in C^), by coordinating their actions and communicating, determine the new feasible state that 
minimizes J(q), solving the optimization subproblem in which all the nodes that are not in C^ t ) keep their 
states constant; 

3) the agents in Cj( t ) actuate the system by updating their state (the injected reactive power). 

In order to describe more precisely this procedure, we rewrite the optimization (9) distinguishing the con- 
trollable and the uncontrollable components of the vector q. Assume with no loss of generality that the first m 
components of q are controllable (i.e. they describe the reactive power injected into the grid by the compensators) 
and that the remaining n — m are not controllable. Let us then partition q as 



q = 



qc 

qv\c 



\ According to this partition of q, let us also partition the matrix Re(X) as 



where q c <G R m and q v \ c E IB 

M N 
N T Q 

Let us also introduce some convenient notation. Consider the subspaces 



Rc(X) 



(10) 



S t :={q C e R m : ]T = , - Vj £ d 

jeCi 



and the m x m matrices 



n 



w 1 



h.keC 



m 



hMeCi 



= diag(l Ci )-— l Ci l£. 

Notice that Si = spanfi^. We list some useful properties of these matrices. 

i) fif = Cli and Vl\ = Qi (where jj means pseudoinverse). 

ii) The matrices Qj, VliMVLi and (fijMOj)" have entries different from zero only in position h, k with 

h,k € Ci. 

iii) ker(Q l Mfi i ) s = kcrfii and span(fi i Mfi i ) tt = spanf^i, and thus (OjMfij)" = (f2jMOj)"f2j. 
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If the system is in the state q and the nodes in the cluster d are activated, they perform the following 
optimization 

argmin J (q' c , q V \c) ■ (H) 
q' c &qc +Si 

Using the standard formulas for quadratic optimization it can be shown that 



q ^ 1 := argmin J(q' c ,q V \c) 



= q c -{Sl i Mn i )\Mq c +Nq vxc ) 
= q c - (fiiMftO'VJ, 



(12) 



where 



(13) 



V J = Mg c + Nq V \ C = [Re(X)q] c e R 

is the gradient of J(qc,qv\c) w i m respect to the decision variables qc- 

We prove now that the computation of q^F 1 ' 1 can be performed by the nodes in a distributed way via a 

suitable approximation. Notice first that, if h £ d, then q^ 1 ' 1 

L J i 

- J2 [(^ Mfi i 



= qh If instead h e Cj, we can write 



opt, 

<7c 



= % 



feec 



From the previous properties we get that 



opt, i 



[n.Mnif] = for all k £ C h and so 

J hk 

= q h -Y] UniMsii)*] [vj] fc . (14) 

In the following, we show how each node h e Cj can compute (14) on the basis of the information available 
to the cluster d. 



B. Hessian estimation from local topology information 

Agents belonging to the cluster d can infer some information about the Hessian M. More precisely they 
can determine the matrix f^MSli appearing in (14). Define now R eS as a m x to matrix with entry in h, k 
equal to Re (Zf^j, where Zf f k are the mutual effective impedances between pairs compensators h, k. From (6) 
and (10) we have that 

Rf k = Rc (Zf k ) = (l h l k ) T M(l h - l fc ) 

and so we can write 

i? eff = - w T M(i h - i fc )ifc 

h,keC 

= diag(M)ll T + 11 T diag(M) - 2M 
where diag(M ) is the to x to diagonal matrix having the same diagonal elements of M. Consequently, since 
fi Ci l = 0, we have that f2ji? eff Oj = -2QjMO^ and then 

n.Mn, = -^n.R^n,. (15) 

We can therefore argue that the elements of the matrix QiMQi can be obtained by the nodes from the mutual 
effective impedances Zf f k for h, k £ d- These impedances are assumed to be known by the nodes belonging 
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to the cluster. This is a reasonable hypothesis, since the mutual effective impedances can be obtained from 
either a priori knowledge of the local microgrid topology, or via online estimation procedures as in [26], or 
via ranging technologies over power line communications as suggested in [27]. 

Then the nodes need to compute the pseudo-inverse of fijMfij. Notice that all these operations have to be 
executed offline once, as these coefficients depend only on the grid topology and impedances. 

C. Gradient estimation via local voltage measurement 

Assume that nodes in d can measure the grid voltage at their point of connection (via phasor measurement 
units that return both amplitude and phase of the measured voltage - see [28], [29]). Let moreover assume the 
following about power line impedances. 

Assumption 3. All power lines in the microgrid have the same inductance/resistance ratio, i.e. 

Z = e j6 Z 

where Z is a diagonal real-valued matrix, whose elements are Z ee = \z e \. Consequently, L = e~i e A T Z~ l A, 
and X := e~i e X is a real-valued matrix. 



This assumption seems quite reasonable in most practical cases (see for example the IEEE standard testbeds 
[30]). The possible effects of this approximation will be commented in Section VII. 
Let then each agent k e d compute 



.(0 



tt ^ \u v \\u k \s'm(Zu k - Zu v 



(16) 



Observe that 



By using (7), we have 



W T 

v k = Im 



-j<> 



-y 

r.A ^ 



U v Uk 



vECi 



(») t 

v).' = Im 



Im 



< "'[ + E>b) {e^l T k Xi + U ) 



1 



~— lc-Xili Xi + e- j2e —- l^XiUo 



\Ci 



+ \ T k XS+\ T k X5 + e^ e \U \ 



where we used the fact that, by (8), U i = s + 5. The third term of the sum is, according to (13), 



Im [l T k Xs] = -l T k Xq = 
The first and the fourth terms of the sum form the quantity 



1 



cos 6 



-(*) T 

v k := Im 



i j 



-,Xi\\Xi + l{X6 



which is infinitesimal for large Uq and so it can be ignored. The second and the last terms form the quantity 



a {i) := Im 



-i*>c± XiUo + e -ie m 



which does not depend on the node h but only on the selected cluster d. Therefore we can write 



(17) 
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D. Description of the algorithm 

From the above consideration, starting from equation (14), we propose the following iterative algorithm. If 
we assume that at time t the system is in the state q(t) £ R™ and that the cluster i is activated. Then the new 
state q(t + 1) will be such that qu{t + 1) = qh(t) for all h £ d, while the node h £ d will inject the new 
reactive power 

Qh(t + 1) = q h (t)- 2 cos 9 Y [(^i? eff ^) J l i/«(t), (18) 

r — * L J hk 



kec t 

where we used that fact that, as ker(fi i i? eff ri i ) l) = kerfii, the term in (17) disappears, while the terms 



have been neglected since they are infinitesimal. 

We have therefore shown as the proposed algorithm can be implemented by the agents of the microgrid in 
a distributed way. In a preliminary, offline phase, each cluster computes (r2ii? eff f2i)". Then, at every iteration 
of the algorithm: 

• a cluster d is randomly chosen; 

• every agent h not belonging to the cluster d holds its injected reactive power constant; 

• every agent h belonging to the cluster d senses the grid voltage at its point of connection, computes 
as in (16), and then updates its injected reactive power according to (18). 

Remark. It is important to notice that the distributed implementation of the algorithm is only possible because 
of the specific physical system that we are considering. By sensing the system locally, nodes can infer some 
global information on the system state (namely, the gradient of the cost function) that otherwise would depend 
on the states of all other agents. For this reason, it is necessary that, after each iteration of the algorithm, 
the agents actuate the system implementing the optimization step, so that subsequent measurement will reflect 
the state change. This approach resembles some other applications of distributed optimization, e.g. the radio 
transmission power control in wireless networks [31 ] and the traffic congestion protocol for wired data networks 
[32], [33]. Also in these cases, the iterative tuning of the decision variables (radio power, transmission rate) 
depends on congestion feedback signals that are function of the entire state of the system. However, these signals 
can be detected locally by each agent by measuring error rates, signal-to-noise ratios, or specific feedback 
signals in the protocol. 

V. Convergence of the algorithm 

In this section we will give the conditions ensuring that the proposed iterative algorithm converges to the 
optimal solution, and we will analyze its convergence rate. 

A. A necessary condition for convergence 
Consider the discrete time system 

q(t + 1) = q{t) - (QiMQif (Mq(t) + Nq v \ c ). (19) 

This system describes the ideal iterative optimization algorithm based on the steps described in (12). It therefore 
approximates the iterative update law (18), and we refer to it for studying the convergence of the proposed 
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algorithm. We introduce the auxiliary variable x = qc—q^ 1 € 1"', where q opt is the solution of the optimization 
problem (9), and q opl are the elements of q opl corresponding to the nodes in C. 

In this notation, it is possible to explicitly express the previous discrete time (19) system as 

x(t + 1) = F l{t) x(t), x(0) e kerl T , (20) 

where 

Fi = I - (QiMQi^M. 

The matrices Fj have some nice properties which are listed below: 

i) the matrices Fi are projection operators, i.e. Ff = F i7 and they are orthogonal projections with respect to the 
inner product (•, -) M , defined as (x, y) M :— x T My. In other words, (Fix, FiX—x) m = x T M(FiX-x) = 0. 

ii) The matrices Fi are self-adjoint matrices with respect to the inner product (•, •) m, i-e- Ff M — MF i7 thus 
having real eigenvalues. 

Notice that, as span(ri l MO i ) tt = spanf^, ker 1 T is a positively invariant set for (20), namely F(kerl T ) C 
ker 1 T for all i. 

It is clear that q{t) converges to the optimal solution q opt is if and only if x(t) converges to zero for any 
initial condition x(0) e ker 1 T . A necessary condition for this to happen is that there are no nonzero equilibria 
in the discrete time system (20), namely that the only x <G ker 1 T such that Fix = x for all i must be x = 0. 
The following proposition provides a convenient characterization of this property. We need the following lemma 
first. 

Lemma 4. x = is the only point in ker 1 T such that Fix = x for all i if and only if 

span[fti . . . ft^] =kerl T . (21) 

Proof: Let us prove the reverse implication first. Assume that span[Oi . . . O^] = ker 1 T and that x G ker 1 T 
is such that Fix = x for all i. Then we can find t/j e R™ such that 

x = ^fliVi- 

i 

Moreover, as Fix = x for all i, then (f^MO^'Mx = and so, by the properties mentioned above, we can 
argue that Mx € kerO^. We therefore have 

x T Mx = vT^Mx = 0, 

i 

and so, since M is positive definite, it yields that x = 0. 

We prove now the other direction. The inclusion span[£li . . . fig] C ker 1 T is always true and follows from 
the fact that spanf^ C kcrl T . We need to prove only the other inclusion. Suppose that x = is the only 
point in ker 1 T such that Fix — x for all i. This means that 

ker 1 T n kcr(7 - Fi) n • • • H ker(7 - F t ) = {0}, 

which implies that 

W n = span 1 + span(7 - Ff) + ■ ■ ■ + span(7 - Ff) 
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Take now any x € ker 1 T . Then there exist aaeR and e R m such that 

e e 



Then al = Mw where 



w = x - y^^fl i Mn i ) i y l =x-'Y^ j CLiZi 



where we have used the fact that span^XSli)" = spanfii. From the fact that spanf^ C kcrl T , we can 
argue that w £ ker 1 T and so it follows that = w T al — w T Mw. Finally, since M is positive definite, we 
can conclude that w = and so that x e span[Oi . . . Q,(\. ■ 
The condition expressed in Lemma 4 can also be expressed as a connectivity test. Let H be the hypergraph 
having the elements in C as nodes and the sets Cj, i = 1, 2, . . . ,£ as hyperedges. This hypergraph describes the 
communication network needed for implementing the algorithm. 

Proposition 5. x = is the only point in ker 1 T such that Fix — x for all i if and only if the hypergraph % 
is connected. 



Proof: Consider the matrix W € K mxm with entries Whk equal to the number of the sets Ci which contain 
both h and k. Let Qw be the weighted graph associated with W. It is quite easy to see that the hypergraph 
having Ci as edges is connected if and only if Qw is a connected graph. 

Let us define 6^ ■ C — > {0, 1} as the characteristic function of the set Ci, namely a function of the nodes 
that is 1 when the node belongs to d and is zero otherwise. Then the Laplacian Lw of Qw can be expressed 
as follows 



l w = ]T (i fc - ifc)(i h - i k ) T w hk 



h.keC 



(ik-ifc)(i(.-i*) T E^(%w 



h,keC 



]T ]T (i h -i fc )(i/.-ifc) T = E 2 l c *l n * 

j=l h,keCi i 



[fii.-.n/] diag{2|Ci|/,...,2|C/|/} 



Notice that the graph connectivity is equivalent to the fact that ker Lw = spanl and so, by the previous 
equality, it is equivalent to 



ker 



Mi 



= span 1 



which is equivalent to (21). The proposition statement then follows from Lemma 4. 
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B. Convergence and rate of convergence 

In order to prove that the connectivity of the hypergraph H is not only a necessary but also a sufficient 
condition for the convergence of the algorithm, we introduce the following assumption on the sequence i(t). 

Assumption 6. The sequence i(t) is a sequence of independently, identically distributed symbols in {1,. ..,£}, 
with non-zero probabilities {pi > 0, i = 1, . . . ,1}. 

Let v(t) := E [x T (t)Mx(t))] = J(q(t)) — J(g opt ) and consider the following performance metric 

R= sup limsup-t^)} 1 /' 

x(0)£kcrl T 

which describes the exponential rate of convergence to zero of v(t). It is clear that R < 1 implies the exponential 
mean square convergence of qc(t) to the optimal solution q^ 1 . 
Using (20), we have 

v(t) = E [x(t) T Mx(t)] = E [x(t) T nMilx(t)] 



'J(0) ' ' ' F o(t-l)^ 



= a 



= .x(0) T E Fj (0) • • • FTQMQF^) ■ ■ ■ F a{0) x(0). 



Let us then define 

A(t) = E [Fj (0) • • • ^.^MnVi) • • • ^(0) • 
Via Assumption 6, we can argue that A(t) satisfies the following linear recursive equation 

A(t+l)=£(A(t)) 

(22) 

A(o) = nMn, 

where £(A) := E [F?AFi\. Moreover we have that 

v(t) = x(0) T A(t)x(0). (23) 
Equations (22) and (23) can be seen as a discrete time linear system with state A(t) and output v(t). 

Remark. The discrete time linear system (22) can be equivalently be described by the equation 

vec(A(t + 1)) = Fvec(A(i)) , 
where vcc(-) stands for the operation of vectorization and where 

F := E [F? ® Ff] e R m2xm2 . 

Notice that F is self-adjoint with respect to the inner product (-, ^m-^m- 1 awc ^ 50 F ami consequently C, 
has real eigenvalues. 

Studying the rate R (which can be proved to be the slowest reachable and observable mode of the system 
(22)-(23) [34]) is in general not simple. It has been done analytically and numerically for some special graph 
structures in [34]. It is convenient to analyze its behavior indirectly, through another parameter which is easier 
to compute and to study. Let us define 

/3 = max{|A| | A e A(F ave ), A ^ 1}, 
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where F ave :=E[F& 

We have the following result. 

Theorem 7. Assume that Assumption 6 holds true and that the hypergraph H is connected. Then 

R<P<1. 



Before proving Theorem 7, observe that the following result descends directly from it. 

Corollary 8. Assume that Assumption 6 holds true and that the hypergraph % is connected. Then the state of 
the iterative algorithm described in Section IV-A converges in mean square to the global optimal solution. 

Remark. Observe that, by standard arguments based on Borel Cantelli lemma, the exponential convergence 
in mean square implies also almost sure convergence to the optimal solution. 

In order to prove Theorem 7 we need a few technical lemmas. 

Lemma 9. Let P,Q € R mxm and P > Q. Then C\P) > C l (Q) for all t e Z> . 

Proof: From the definition of £, we have 

x T [C(P) - C(Q)\ x = x T [E [FfPFi] - E [F^QF,]] x 

= E [x T FT (P-Q)F iX \ > 0. 
By iterating these steps t times we then obtain £'(-P) > C l {Q). ■ 

Lemma 10. For all A we have that IXC'^Afi)^ = JXC*(A)ft. 

Proof: Proof is by induction. The statement is true for t = 0, as il 2 = ft. Suppose it is true up to t. We 
then have 

n£ t+1 (A)n = Q£(£*(A))ft = fi£(fi£*(A)fi)0 

= ncinc^nAn^n = nc t+1 (riAn)n. 



Lemma 11. All the eigenvalues of F ave are real and have absolute value not larger than 1. If span fl 1 - 
ker 1 T and if Assumption 6 holds, then the only eigenvalue of F ave on the unitary circle is A = 1, with multiplicity 
1 and with associated left eigenvector 1 and right eigenvector M _1 1. 

Proof: Recall that the matrices Fi are projection operators, i.e. Ff = F{ and so they have eigenvalues 
or 1. Recall moreover that Fi are self-adjoint matrices with respect to the inner product (•, -)m, defined as 
(x, j/)m = x 7 My, This implies that H-F^Hm < 1 where || • ||m is the induced matrix norm with respect to the 

1/2 

vector norm ||x||m := (x, x )m ■ This implies that 

II -Fave || M = 

||E[F 4 ]|| M <E[||^|| M ] < 1 



October 24, 201 1 



DRAFT 



18 



And so, since also F ave is self-adjoint, its eigenvalues are real and are smaller than or equal to 1 in absolute 
value. 

Assume now that = Ax, with |A| = 1. Then we have 

\\x\\ M = ||Fave||M = ||E [Fi] x\\ M < E [\\Fix\\ M ] < \\x\\m- 

If Assumption 6 holds (and thus the probabilities p^'s are all strictly greater than 0), then the last inequality 
implies that H-F^xHm = \\ x \\m for all i's. As Fi are projection matrices, it means that Fix = x and so 
Mx e ker flf , Vi. Using the fact that span fl 1 ■ ■ -fl e — ker 1 T , we necessarily have x = M _1 1. By inspection 
we can verify that the left eigenvector corresponding to the same eigenvalue is 1. ■ 
We can now give the proof of Theorem 7. 

Proof of Theorem 7: Let us first prove that f2£(f2Mfi)fi < /3QMQ. Indeed, we have 

x T VLC{VtMQ,)Vtx = E [x T flFTnMflFiQx\ 
= E [x T 9.Fj'MF i 9.x\ 

= x T nM 1 / 2 E [m^f.m- 1 ' 2 ] M^Ctx, 

where we use the fact that JlFjfi = F t fl and Fj MFi = MF t . Notice moreover that E [M 1 / 2 F l M~ 1 / 2 ] = 
M x l 2 F W zM~ x l 2 is symmetric and, by Lemma 11, it has only one eigenvalue on the unit circle (precisely in 
1), with eigenvector M~ 1/2 1. As M 1/2 ^lx _L M" 1/2 1 for all x, we have 

x T Q.C{Q.MQ)Q.x < /3x T nMflx, 

with (3 = max{|A| | A e A(F ave ), A ^ 1}. 

From this result, using Lemmas 9 and 10, we can conclude 

Q£*(OMfi)n = SIC*- 1 (£(OMfi)) n 

= nc 1 - 1 (0£(OMn)fi) n 

= pric*- 1 (omo) n < ■ ■ ■ < frriMn, 

and therefore R< (3. ■ 
The tightness of (3 as a bound for R has been studied in [34], by evaluating both j3 and R analytically and 

numerically for some special graph topologies. In the following we consider (3 as a reliable metric for the 

evaluation of the algorithm performances. 

The following result shows what is the best performance (according to the bound j3 on the convergence rate 

R) that the proposed algorithm can achieve. 

Theorem 12. Consider the algorithm (19), and assume that % describing the clusters Ci is an arbitrary 
connected hypergraph defined over the nodes C. Let Assumption 6 hold. Then 

(Ef=ift|Ci|)-i 

/?>1- A - 1 • (24) 

m — 1 

In case all the sets Ci have the same cardinality c, namely \d\ = c for all i, then 

P > 1 - — \- (25) 

m — 1 
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Proof: Let 

E t := (QiMfl^M (26) 

e 

E ave :=E[E i ]=J2piE i (27) 
»=i 

It is clear that /? = 1 - /?', where /?' = min{|A| | A e A(_E ave ), A 7^ 0}. We have that 

J2 A, = Tr(£ ave ) = Tr ( £ j = £ Tr . 

Aj-eA(£; av e) \«=i / «=i 

Notice now that i?i have eigenvalues zero or one and so 

Tr (Ei) = rank Ei = dim (span Ei) = dim (span(fi i Mfi i ) ti ) 
= dim (span f^) = |Ci| — 1. 



This implies that 



and so 



E x i = (E^i) - 1 ' 

Aj-eA(£; ave ) \i=l / 

£' = min{|A| | Ae A(£ ave ),A^0} 

< EA jg A( g „ vc )A 3 = (ELiPd^) -1 

m — 1 m — 1 



VI. Optimal communication hypergraph for a radial distribution network 

In this section we present a special case in which the optimal convergence rate of Theorem 12 is indeed 
achieved via a specific choice of the clusters C,. To obtain this result we need to introduce the following 
assumption, which is commonly verified in many practical cases (see for example the standard IEEE testbeds 
[30]). 

Assumption 13. The distribution network is radial, i.e. the corresponding graph Q is a tree. 
We start by providing some useful properties of the matrices Ei. 

Lemma 14. The following properties hold true: 

i) span Ei = span f^; 

ii) Ef = E t ; 

Hi) EiX — x for all x <E spanfij. 
Proof: 

i) First observe that, since M is invertible, then span E t = span(fi l M^ i )»M = span(O i M^)«. Finally 
from the properties of the matrices (f^MQj)", we obtain that spani^ — spanf^. 



October 24, 201 1 



DRAFT 



20 



ii) First observe that E? = (I - F,) 2 = I - + Ff. Now, since Ff = F t we obtain that I - 2F, + Ff = 
I-F i = Ei. 

iii) Let x G spanOi = span .5$. Then x = EiZ for some vector z and so EiX = Efz = EiZ = x. 

■ 

Let us now consider for any pair of nodes h, k the shortest path V hk C £ (where, we recall, a path is a 
sequence of consecutive edges) connecting h, k. Define moreover 

v Ci ■■= IJ v hk . 

Lemma 15. Let Assumption 13 hold true. Then 

i) (l fc - l fe ) T M(l r - 1.) = Eeev hk nv rs M^)- 

ii) Ei(l h - l k )=0 if v Ci nv hk = <D. 

Proof: 

i) Observe that (1^ — l k ) T M(l r — l s ) = (1/, — lfc) T X(l r — l s ), where with some abuse of notation we 
have denoted with the same symbol 1/, the vector in K m and the corresponding vector in W 1 . Notice that, 
if we have a current i = l r — l s in the network, then according to (7) we get a voltage vector 

u = Xi + 1U = X(l r - l s ) + 1U 

and so (1/j — lfe) T X(l r — l s ) = Uh — u k . Observing the following figure 




h S ^ k 

we can argue that u h -u k = u h > - u k > = J2eev h , k , z e and so - lfc) T M(l r - l s ) = J2eev h , k , R- e ( z e)- 
ii) If V Cz n V hk = 0, then 



fW(ih-ifc) 



2\C 



l — (!r - l-)(lr - ls) T M(U - lfc) = 



since V hk nT rs = ® for all r, s e Cj. Thus £;(!,, - l fe ) = (O^O^O^l,, - l fc ) = 0. 



We give now the key definition to characterize, together with connectivity, the optimality of the communication 
hypergraph H. 

Definition 16. The hyperedges {Ci, i = 1, 2, . . . ,£} of the hypergraph H are edge-disjoint 1 if for any i,j such 
that i 7^ j, we have that 

V Ci n V Cl = 0. 
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In order to prove the optimality of having edge-disjoint clusters, we need the following lemma. 

Lemma 17. Assume that Assumptions 13 holds true and that the clusters Ci are edge-disjoint. Then EiEj = 
for all i 7^ j. 

Proof: Let x be any vector in R m . Notice that, since sp&nEj = spanfij, then there exists z S R m such 

that 

EiEjX = EiQjz = ^— E t (l r -l s )(l r -l s fz = 

where the last equality follows from the fact that Vci H V rs = for all r,s € Cj. ■ 
From the previous lemma we can argue that every vector in spanfij is eigenvector of E. dwe = ^f=i Pi^i 
associated with the eigenvalue pi. On the other hand, recall that the hypergraph connectivity implies that 
kerl T = Y^i=i spanili = Ym=i spanEj and hence we have that, besides the zero eigenvalue, the eigenvalues 
of E. dve are exactly p\,...,pt and so 

P = 1 - min { Pi \. 

The bound (3 is minimized by taking pi = \ ji for all i, obtaining 

" _1 7 

Notice finally that E^Ej = for all i ^ j implies that J2i=i spanEi is a direct sum and then 



m 



1 = ^ dim span ^ =^(|C,| - \)= J £ i \C i \-L 



If all the sets d have the same cardinality c, then Ic = m + 1 — 1, and so 



I m — 1 

which, according to Theorem 12, shows the optimality of the hypergraph. 

VII. Simulations 

In this section we present numerical simulations to validate both the model presented in Section III, and the 
randomized algorithm proposed in Section IV. 

To do so, we considered a 4.8 kV testbed inspired from the standard testbed IEEE 37 [30], which is an actual 
portion of power distribution network located in California. We however assumed that load are balanced, and 
therefore all currents and voltages can be described in a single-phase phasorial notation. 

As shown in Figure 2, some of the nodes are microgenerators and are capable of injecting a commanded 
amount of reactive power. Notice that the PCC (node 0) belongs to the set of compensators. This means that 
the microgrid is allowed to change the amount of reactive power gathered from the transmission grid, if this 
reduces the power distribution losses. The nodes which are not compensators, are a blend of constant-power, 
constant-current, and constant-impedance loads, with a total power demand of almost 2 MW of active power 
and 1 MVAR of reactive power (see [30] for the testbed data). 

'This definition of edge-disjoint for this specific scenario resembles a similar concept for routing problems in data networks, as for 
example [35]. 
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Figure 2. Schematic representation of the IEEE37 testbed. Circled nodes represent microgenerators taking part to the distributed reactive 
power compensation. The thick curved lines represent clusters of size c = 2 (i.e. edges of H) connecting pair of compensators that can 
communicate and coordinate their behavior. 

The length of the power lines range from a minimum of 25 meters to a maximum of almost 600 meters. The 
impedance of the power lines differs from edge to edge (for example, resistance ranges from 0.182 O/km to 
1.305 51/km). However, the inductance/resistance ratio exhibits a smaller variation, ranging from Zz(e) = 0.47 
to Zz(e) = 0.59. This justifies Assumption 3, in which we claimed that Zz(e) can be considered constant 
across the network. The effects of this approximation will be commented later. 

Without any distributed reactive power compensation, distribution power losses amount to 61.6 kW, 3.11% 
of the delivered active power. 

Given this setup, we first estimate the quality of the linear approximated model proposed in Section III. As 
shown in Figure 3, the approximation error results to be negligible. 

We then simulated the behavior of the algorithm proposed in Section IV. We considered the two following 
clustering choices: 

• edge-disjoint gossip: motivated by the result stated in Section VI, we enabled pairwise communication 
between compensator in a way that guarantees that the clusters (pairs) are edge-disjoint; the resulting 
hypergraph % is represented as a dashed line in Figure 2; 

• star topology: clusters are in the form = {0,v} for all v e C\{0}; the reason of this choice is that, 
as is the PCC, the constraint l T q = is inherently satisfied: whatever variation in the injected reactive 
power is applied by v, it is automatically compensated by a variation in the demand of reactive power 
from the transmission grid via the PCC. 

In Figure 4 we plotted the result of a single execution of the algorithm in the edge-disjoint gossip case. One 
can see that the algorithm converges quite fast, reducing losses to a minimum that is extremely close to the 
best achievable solution. The results achieved by the proposed algorithm on this testbed are summarized in the 
following table. 
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Figure 3. Comparison between the network state (node voltages) computed via the exact model induced by (1), (2), (4), and (3) (o), and 
the approximate model proposed in Proposition 2 (+). 




Figure 4. Power distribution losses resulting by the execution of the proposed algorithm. A edge-disjoint hypergraph, yielding optimal 
convergence speed, has been adopted. The dashed line represent the minimum losses that can be achieved by properly choosing the reactive 
power injected by the compensators. 



Losses before optimization 61589 W 

Fraction of delivered power 3.11 % 

Losses after optimization 50338 W 

Fraction of delivered power 2.55 % 

losses reduction 18.27 % 

Minimum losses J opt 50253 W 

Fraction of delivered power 2.54 % 

losses reduction 18.41 % 

The minimum losses J opt has been obtained via numerical optimization performed on the exact microgrid 
model, and represent the minimum losses that can be achieved by properly choosing the amount of reactive 
power injected by the compensators (and retrieved from the PCC). The difference between this minimum and the 
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Figure 5. Simulation of the expected behavior of the algorithm for different clustering strategies. The difference between power losses 
at each iteration of the algorithm (averaged over 1000 realizations) and the minimum possible losses, has been plotted for two different 
clustering choices: edge-disjoint gossip (solid line) and star topology (dashed). The dotted line represent the best possible performance. 



minimum achieved by the algorithm proposed in this paper is partly due to the approximation that we introduced 
when we modeled the microgrid (assuming large nominal voltage Uq), and partly due to the assumption that 
6 is constant across the network. Further simulative investigation on this testbed showed that the effect of this 
last assumption is largely predominant, compared to the effects of the approximation in the microgrid state 
equations. Still, the combined effect of these two terms is minimal, in practice. 

In Figure 5 we compared the behavior of the two considered clustering strategies, plotted together with the 
curve representing the best achievable performance as given in Theorem 12. One can notice different things. 
First, both strategies converge to the same minimum, which is slightly larger than the minimum losses that 
could be achieved by solving the original, nonconvex optimization problem. As expected, the clustering choice 
do not affect the steady state performance of the algorithm. Second, the performance of the edge-disjoint gossip 
algorithm results indeed to be better than the star topology in the long-time regime, as the analysis in Section VI 
suggests, and corresponds to the fastest achievable performance, as predicted by Theorem 12. Third, one can see 
that, while the performance metric that we adopted is meaningful for describing the long-time regime, it does 
not describe the initial stage (or short-time behavior), when the full dynamics of the system (19) contribute to 
the algorithm behavior. This fact suggests that a different metric should possibly be adopted, in order to better 
characterize the two regimes. The choice of the most appropriate metric requires a preliminary characterization 
of how fast the reactive power demands change in the microgrid, compared with the rate of execution of the 
algorithm iterations. If the time-varying demands are much slower than the algorithm, then we expect that the 
analysis of the slowest eigenvalue remains the most meaningful choice. On the contrary, if the rate of change 
of demands is comparable with the rate of execution of the algorithm, then we expect that a different index 
will have to be considered to characterize the algorithm performance. 

VIII. Conclusions 

The proposed model for the problem of optimal reactive power compensation in smart microgrids exhibits 
two main features. First, it can be casted into the framework of quadratic optimization, for which robust solvers 
are available and the performance analysis becomes tractable; second, it shows how the physics of the system 
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can be exploited to design a distributed algorithm for the problem. 

On the basis of the proposed model, we have been able to design a distributed, leaderless and randomized 
algorithm for the solution of the optimization problem, requiring only local communication and local knowledge 
of the network topology. 

We then proposed a metric for the performance of the algorithm, for which we are able to provide a bound 
on the best achievable performances. We are also able to tell which clustering choice is capable of giving 
the optimal performances. It is interesting that the optimal strategy requires short-range communications: this 
is somehow surprising, considered that consensus algorithm (which share many features with the proposed 
algorithm) benefit from long-range communication that shorten the graph diameter. At the same time, this 
result is also motivating from the technological point of view, because achieving cooperation of inverters that 
are far apart in the microgrid presents a number of issues (e.g. time synchronization, limited range of power 
line communication technologies). 

Motivated by the final remarks of Section VII, we plan to add two elements to the picture in the next future 
in order to understand if the proposed metric is really meaningful for the real problem. These elements are the 
network dynamics (both the dynamic behavior of the power demands, and the electrical response of the system 
when actuated), and the presence of operational constraints for the compensators, whose capability of injecting 
reactive power is limited and possibly time varying. 

Appendix 
Proof of Proposition 2 

We first introduce the new variable e := 1/Uq. This way we can see the currents i and the voltages u as 
functions i(e),u(e) of e, determined by the equations 

u = Xi + e- 1 ! 

< l T i = (28) 

u v i v = s v \eu v \ Vv , Vu e V\{0} 
We want to give the Taylor approximation of i(e) around e = 0. 
Let us define 

5 v (e) := i v (e)e _1 - s v 
X v (e) := u v (e)e - 1 

and substitute them into (28). We get a system of equations in 6, A, e which can be written as 

G v (6,X,e)=0, VueV 

< 

F v (6,X,e) = 0, VveV 

where 

Go(S, A, e) := A 

G V (S, A, e) := A.„ - J2 we v X vw(s w + <5„,)|e| 2 , Vw ^ 

(29) 

F (5,\,e) := J2 we v S *> 

F V (S, A, e) := (1 + \ v )(s v + 5 V ) -s v \l + A„|"" , Vv ± 
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We have 2n complex equations in 2n + 1 complex variables. 

We interpret now the complex numbers as vectors in R 2 and G v , F v as functions from 



to 



Observe that G(S, A, e)\s=o,\=o, e=o = and A, e)\s=o,X=o,e=o = for all v e V. By the implicit 

function theorem, in order to prove that 6(e), A(e) are infinitesimal in e, it is enough to prove that the matrix 

'dG v \ (dG v s 



d5 u 



w / v.wev 



dF v 
86, 



d\ w / ,. .,. 
( dF v 

\dX w J v weV 



evaluated in 5 = 0, A = 0, e = is an invertible. 

In order to determine this matrix we need to introduce the following notation. If a £ C, and a = a R + ja 1 
with a 7 *-, a 7 G R, then we define 



(a) := 



a fl -a 7 



a J a 

With this notation observe that, if a, x e C and we consider those complex numbers as vectors in R 2 and if 
we consider the function j(x) := ax as a function from R 2 to R 2 , then we have that df/dx = (a). Notice 
moreover that the the function g(x) := x has 

r 1 

-1 



8g 
dx 



N :-- 



From these observations we can argue that 

dGo 



d5 u 

dFo 
d5 w 



= 0, 



dG 
d\ w 

dFp 
d\ w 




Instead for all v ^ we have that 



while 



dG v 
dX m 




dGv 
d5 v , 



= (-X^|e| 2 ) VweV 



which evaluated in e = yields 0. Observe finally that, if w ^ v, then dF v /dX w = dF v /d5 w = and that 



dF, 
d5 v 



v - = (1 + X V )N 



which evaluated in = yields the matrix TV. On the other hand, 

dFv_, d\l + X v \^ 

dx v ~ {Sv + v) Sv dx v 

where s v is a 2-dimensional column vector and is a 2-dimensional row vector. With some easy 



computations it can be seen that 

d\l + X v \ v " 
dX v 



n v [(i + Af f + (A 7 ) 



>7\21 — 



1 + Xy A 7 
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which evaluated in A^ = yields [rj v 0]. Hence in A„ = we have 

dFy 



Qv ■= (s v ) - Sv[r)v 0] 



We can conclude that for S = 0, A = 0, e = we have 



(dG v 



dF v 
86,., 



v,wEV 



which is invertible. 



dG v 



\ 96 w J v weV \d^wJv,wev 
dF v 



<9A„ 



w / v,wEV 






• 


• 


I 


• 








• 


• 





I ■ 








• 


• 





• 


I 


/ 


I ■ 


• I 





• 








N ■ 


■ 





Qi ■ 








• 


• N 





• 


• Qn-1 
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